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Pref ace 

This  report  contains  material  that  was  distributed 
at  the  Workshop  on  ID  and  2D  Transport  Codes  which 
met  at  Cornell  University  in  August  1978.   This  is 
a  companion  to  Report  MF-91,  "Classical  Diffusion: 
Theory  and  Simulation  Codes",  March  1978,  which 
describes  the  analytic  background  for  the  codes 
described  here.   The  persons  responsible  for  creation 
of  the  individual  codes  are  named  in  the  cited  ref- 
erences and  are  primarily  D.  C.  Stevens,  E.  Turkel, 
and  P.  N.  Hu  with  the  assistance,  in  later  codes, 
of  W.  Grossmann  and  S.  Wollman. 
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Abstract 

A  survey  is  given  of  a  family  of  classical  transport 
codes,  recently  termed  "lyD",  which  efficiently  and 
accurately  follow  the  evolution  of  plasma  configura- 
tions on  a  long  time  scale,  following  coupled  changes 
in  plasma  shape  and  topology  with  transport  (but  not 
wave  motion) .   Codes  have  been  constructed  and 
operated  (since  1974)  which  include  various  combina- 
tions of  finite  beta,  general  plasma  cross-section 
and  aspect,  various  topologies  (Doublet,  tearing, 
reversed-f ield  mirror)  including  time  dependent 
transitions  in  topology  resulting  from  external  coil 
variation  and  plasma  transport,  with  models  including 
(classical)  tensor  resistivity  and  heat  flow  as  well 
as  the  adiabatic  limiting  case. 
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Introduction 

The  class  of  codes  which  have  recently  been  termed  "lyD" 
are  actually  two  dimensional;  they  are  codes  which  accurately 
solve  generically  2D  (e.g.  axial  and  helically  symmetric)  classi- 
cal transport  (and  adiabatic)  problems .   The  algorithms  involve 
carefully  selected  interactions  among  some  operations  which  are 
ID  and  others  which  are  2D.  The  procedures  were  originally  intro- 
duced (1970)  to  follow  the  time  evolution  of  the  Pf irsch-Schluter 
model  and  were  then  extended  and  developed  to  efficiently  march 
in  the  Grad-Hogan  transport  formulation  in  which  the  time  scales 
of  changes  in  plasma  shape  and  topology  are  comparable  with  and 
coupled  to  the  various  transport  times,  but  in  which  all  this 
activity  takes  place  slowly  compared  to  the  Alfven  transit  time. 

The  characteristic  feature  of  the  lyD  family  of  algorithms 
is  the  separation  of  the  solution  process  into  alternate  marching 
in  ID  of  profiles  (current,  temperature,  flux,  volume,  etc.)  and 
2D  evaluation  of  contours  and  geometrical  coefficients  (surface 
inductance,  surface  averaged  resistivity,  etc.).   The  basic 
present  limitations  of  the  lyD  technique  are  to  macroscopic, 
scalar  pressure,  MHD  models.  Codes  have  been  constructed  and 
operated  (since  1974)  which  include  various  combinations  of  finite 
6,  general  plasma  cross-section  and  aspect,  various  topologies 
(Doublet,  tearing,  reversed-f ield  mirror)  including  time  dependent 
transitions  in  topology  resulting  from  external  coil  variation 
and  plasma  transport,  with  models  including  (classical)  tensor 
resistivity  and  heat  flow  as  well  as  the  adiabatic  limiting  case. 
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The  three  irost  important  advantages  of  these  methods  as 
compared  to  fully  classical  (including  inertia  and  wave  motion) 
and  compared  to  the  very  recently  developed  2D  marching  of  the 
Grad-Hogan  model  (by  S.  C.  Jardin)  are: 

1)  complex  plasma  geometry  and  topology,  including 
dynamic  changes  in  topology  are  easily  managed; 

2)  geometry  changes  and  profile  changes  each  march  with 
its  own  size  time  step; 

3)  the  code  logic  is  such  as  to  be  easily  adaptable  to 
convert  any  ID  (cylindrical)  "Kitchen  Sink"  transport  code  to 
a  fully  2D  basis. 

Note:   conventional  terminology  terms  ID  a  code  in  which 
all  the  various  physical  processes  are  assembled  in  ID  even 
though  individual  processes,  such  as  injection  or  a-particles, 
may  be  treated  in  2D  or  even  3D  (but  not  self-consistently) . 

It  is  expected  that  extension  of  IjD  techniques  to  closed 
line  3D  (EBT)  and  to  completely  general  3D  toroidal  transport 
should  not  present  serious  difficulties.   On  the  other  hand, 
strongly  anisotropic  problems  (such  as  guiding  center  or  drift 
models  for  a  mirror  machine  with  loss  cone)  do  not  seem 
amenable  at  present.   The  neoclassical  (guiding  center  or  drift 
but  approximately  scalar  pressure)  problem  seems  to  have 
reached  a  stage  close  to  that  of  the  classical  MHD  Grad-Hogan 
model  some  five  years  ago;  the  theory  including  consistent 
geometry  changes  seems  to  be  visible;  practical  algorithms  and 
operating  codes  for  a  lyD  or  2D  version  are  still  to  be  implemented. 
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The  following  section  contains  a  brief  summary  of  the 
generic  l^D  algorithm  (not  specialized  to  any  particular 
configuration) .   This  is  followed  by  a  more  detailed  description 
of  six  specific  codes  (each  with  a  few  variations)  and  each 
with  results  from  a  few  typical  runs.   Except  for  the  most 
recent  code,  much  of  the  information  has  been  published,  but 
nowhere  conveniently  assembled  in  one  place.   Our  intent  in  this 
report  is  to  show  the  interconnections  and  common  features  of 
a  large  number  of  plasma  configurations  and  models,  frequently 
accompanied  by  spontaneous  changes  of  topology,  handled  by 
similar  numerical  methods,  but  only  in  enough  detail  to  paint 
a  broad  picture  of  the  methods  and  results. 
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The  Basic  IjD   Algorithm 


One  part  of  the  code  is  an  equilibrium  solver.  In  2D,  given 
p(^)    and  B   =  f(^)  profiles,  the  problem  is  a  nonlinear  elliptic  PDE, 

(1)  Aijj  =  F{^)     ,  F  =  -  dp/dij;  -  fdf/dif; 

One  method  of  solution  is  to  iterate  using  a  Poisson  solver, 

(2)  Aip^_^^  =  F[i|;^(x,y)]  =  ^^i^'V) 

This  may  or  may  not  converge  in  cases  where  the  solution  is 
unique  as  well  as  (after  bifurcation)  where  the  solution  is  non- 
unique.   There  is  a  large  literature  and  a  stock  of  more  and  more 

* 

sophisticated  procedures,   but  there  is  no  such  thing  as  a 

universally  convergent  algorithm.   In  a  case  of  multiple  solutions, 
numerical  stability  of  the  algorithm  (distinct  from  physical 
stability)  determines  which  of  the  solutions  may  be  found. 
A  second  formulation  of  the  equilibrium  problem  is 

(3)  A^l)    =  F(V) 

where  the  RHS  is  a  known  function  of  the  volume,  V,  within  a 
ijj-contour.   The  same  iteration 

(4)  ^'^„  +  i  =  F(V^)  ,  V   =  V  (x,y)  ^  i)    (x,y)  =  const. 

n+i      n     n    n         n 

is  usually  more  rapidly  convergent  than  (1)  (but  requires  contours 
of  ii   and  volumes  to  be  evaluated)  .  The  numerical  stability  of 
solution  algorithms  for  (3)  is  different  from  (1)  and  both  are 
different  from  physical  stability. 

*  See  for  example,  B.  Harder  and  H.  Weitzner,  Plasma  Phys.  12, 

435,  1970,  and  R.  Meyer-Spasche,  Max  Planck  Institute  of  PTasma- 
physics  Report,  IPP  6/141,  December,  1975. 
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A  third  equilibrium  formulation,  given  adiabatic  profiles 
U(^)    and  \)  {4))    defined  by 

Y 


(5) 


P  =  y [^' (V)] 
f  =  vii'  (V) 


leads  (for  y   =    2    and  v  =  0)  to  the  QDE  (or  GDE)   [Ref.  (1)] 


(6) 


A^p    =    -    ^    i^P')^    -    2y4;"  =    R 


In  a  configuration  with  simple 
contours  as  shown,  equilibrium 
equations  (1)  and  (3)  are  well- 
posed  with  a  given  elliptic 
boundary  condition,  iJj  ;  but  the 
QDE  (6)  also  demands  specifica- 


tion of  (ii  ,  at  the  "center", 
^o 


The  QDE  will  amost  never  converge  on  simple  iteration 


(7) 


All;  _^,  =  R 
n+1     n 


(and  this  gives  no  opportunity  to  impose  proper  boundary 
conditions) .   A  method  which  converges  extremely  well  is  to 
alternate  between  the  Poisson  equation 


(: 


All;  ,  ,  =  R  (x,y) 
n+1     n    -^ 


which  is  used  to  compute  only  the  geometric  contours  11^,-1=  const, 
(the   i|;   ,  profile  is  discarded)  and  the  averaged  QDE.   Given  the 
contours,  the  inductance  profile 


-8- 


(.9)  ^n^^^  "  <|  VV|^>  -  I  I  Vv|ds 

is  calculated  and  the  average  pressure  balance  (average  QDE) 

(10)  (Ktl>')  •  ^  -  ^  (i^')^  -  2M'J^" 

is  solved  as  an  ODE  for  ^  {V)  ,    given  y(i|^)  [a  priori]  and  given 
K(V)  [temporarily,  by  iteration].   From  ^ (V)    and  the  previous 
contours,  V(x,y),  one  evaluates  R  ^■,(x,y)  and  repeats  the  pro- 
cedure (which  usually  is  found  to  converge  very  rapidly) . 

The  ODE  allows  specification  of  <|j  at  the  center  (V==0)  as 
well  as  the  elliptic  boundary  condition  at  the  plasma  edge;  after 
convergence  of  the  iterations,  both  ODE  and  PDE  boundary  condi- 
tions are  satisfied;  (the  combination  specifies  a  well-posed 
problem  for  the  QDE)  . 

In  a  nonadiabatic  problem  with  transport,  the  averaged  (ID) 
transport  equations  determine  the  time  variation  of  lJ(i|;)  and 
v(i|j).   In  2D,  a  sample  transport  equation 

(11)  1^  +  u-v^  =  nAip 

o  t 
becomes,  when  averaged, 

(12)  ^^  +   ui|^'  =  n(K4j')  • 
where 

(13)  §1  =  1^  'l^(x,y,t)   ,  ,J;^  =  1^  ^(V,t) 

t 
U  =  <u'VV>  =  ()  u'dS 

If  an  adiabatic  variable  (e.g.  mass)  is  taken  as  the  independent 
variable. 
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ijj^  +  UijJ'  =  di|;/dt  =  1^  ^(M,t)  ,   M  = 


pdV 


then  the  velocity  U  (the  net  plasma  diffusion  velocity)  can  be 
eliminated  and  the  entire  set  of  transport  equations  become 
ID  transport  equations  in  M  and  t,  with  various  geometrical 
coefficients  similar  to  K.   Note  that  U,  which  is  the  goal  of 
classical  diffusion  transport  theory  is  here  eliminated;  (one 
can  go  back,  afterwards,  and  compute  U  =  (d/dt) V (M, t) ) . 

In  a  Pf irsch-Schliiter  formulation  where  U  is  given  by  a  trans- 
port relation,  this  transport  formula  can  be  used  to  eliminate  U 
from  all  equations,  and  the  volume  V  (for  example)  can  be  used 
as  the  independent  ID  variable.   In  the  more  general  Grad- 
Hogan  model,  U  (mass  flow)  is  not  a  transport  variable,  and 
the  ID  formulation  (with  U  eliminated)  is  obtained  only  by  use 
of  appropriate  independent  and  dependent  variables  (usually, 
but  not  always,  adiabatic  variables)  [Ref .  4]. 

There  is  always  one  fewer  averaged  transport  equation 
than  originally  (since  one  of  the  "diffusing"  variables  is 
taken  as  the  independent  ID  variable) .   The  count  is  made  up 
by  adding  the  average  pressure  balance  (10)  [in  an  appropriate 
independent  variable,  not  necessarily  V]. 

The  skeleton  algorithm  is: 
1)    From  initial  equilibrium  calculate  (as  contour  averages) 
coefficients  for  ID  averaged  transport  equations.   These  trans- 
port equations,  together  with  the  average  pressure  balance,  are 
marched  in  ID  for  many  diffusion  time  steps  (to  advance  one 
geometry  time  step) . 
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2)    Using  the  diffused  ID  profiles  at  the  advanced  geometry 
time  step,  a  RHS  R(x,y)  is  computed  for  the  Poisson  equation  to 
update  the  geometry;  [here  one  has  a  choice,  involving  simpli- 
city, accuracy,  stability,  etc.,  of  using  the  forms  (1)  or  (3) 
or  (6)  since  all  profiles  are  available  from  the  result  of 
diffusion] . 

There  are  two  iteration  loops,  the  inner  one  to  solve  for 
a  geometry  given  F{i(j)  or  F  (V)  or  \i(\p)  ;    and  the  outer  one  which 
repeats  the  diffusion  over  a  geometry  time  step  until  all 
geometrical  profiles  such  as  K(V)  converge  at  the  advance  time. 

One  can  be  naive  and  take  K(V,t)  piecewise  constant  in  t 
(no  iteration  of  the  outer  loop) ,  or  make  it  linear  in  t  between 
geometry  steps  (for  insertion  in  the  transport  equations),  or 
more  sophisticated  (e.g.  predictor-corrector).   In  axial  symmetry, 

K  has  to  be  treated  more  carefully  than  the  other  geometrical 

2        2 
profiles,  <l/r  >  or  <r  >.   An  important  point  in  marching  the 

geometry  with  a  separatrix  is  to  write  profiles  in  normalized 

coordinates,  V/V    inside  and  (V-V    ) /V  -V    )  outside.  The 
'  '    Sep  sep"  p   sep 

outer  loop  will,  in  this  case,  also  advance  V    (t) ;  this  is 
quite  sensitive. 

There  are  at  least  three  (and  preferably  four)  independent 
meshes 

1)  2D   -    (x,y) 

2)  ID   -   contours    in    (x,y)    to   evaluate   averages    [coarse] 

3)  ID   -   diffusion    [fine] 

4)  ID   -    interpolation   among   other  meshes    [fine] . 
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There  are  also  two  (partly)  independent  convergence  criteria 
e.  and  e   for  the  inner  and  outer  loops.   The  optimization 
problem  is  formidable  and  depends  on  the  goal.   For  example,  in 
studying  the  small  resistivity  approach  to  adiabatic,  there  is 
a  very  sharp  moving  singularity  in  K(V,t)  at  a  separatrix;  this 
turns  out  to  require  careful  (nonuniform)  contour  placement. 
To  optimize  computer  time,  one  will  play  with  more  sophisticated 
extrapolations  of  K(V,t)  in  t,  with  the  values  of  £.  and  e  , 
with  the  relative  sizes  of  2D,  contour,  and  diffusion  meshes, 
as  well  as  contour  accuracy,  special  handling  near  islation 
and  at  the  plasma  edge,  and  ID  transport  algorithms. 

In  one  slightly  optimized  transition  from  a  Belt  Pinch  to 
a  fully  formed  Doublet  with  practical  transport  parameters,  it 
was  found  to  be  sufficient  to  recompute  the  goemetry  only  five 
times.   In  a  transition  from  circular  to  D-shaped  (or  shift 
induced  by  increasing  3)  it  should  rarely  be  necessary  to  use 
more  than  two  or  three  geometry  calculations;  but  this  will  have 
to  be  confirmed  in  each  given  problem. 

More  details  of  the  equations  and  algorithms  are  found  in 
the  following  publications  [for  additional  references,  see  (4)]: 

1.  Grad,  H. ,  Hu,  P.N.,  and  Stevens,  D.C.  "Adiabatic  Evolution 
of  Plasma  Equilibrium,"  Proc.Nat. Acad. Sci. ,  USA  72 , 

pp.  3789-3793  (1975) . 

2.  Grad,  H.,  Hu,  P.N.,  Stevens,  D.C.,  and  Turkel,  E. ,  "Diffusive 
Transition  from  Belt  Pinch  to  Doublet,"  2nd  European  Conf.  on 
Computational  Physics,  Garching,   Germany,  April  27-30, 
1976. 

3.  Grad,  H.,  Hu,  P.N. ,  Stevens,  D.C.,  and  Turkel,  E.,  "Classical 
Plasma  Diffusion,"  IAEA  6th  Conf.  on  Plasma  Physics  and  CTR, 
Berchtesgaden,  Germany,  Oct.  6-13,  1976. 
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4.  Grad,  H.,  Hu,  P.N.,  "Classical  Diffusion:  Theory  and 
Simulation  Codes",  Workshop  on  High  Beta  Plasma,  Varenna, 
Italy,  Sept.  1977.   Also  Report  MF-91,  CIMS,  NYU,  March 
1978. 

5.  D.  Nelson,  "Intense  Neutral  Beam  Heating  in  the  Adiabatic 
Approximation",  Workshop  on  High  Beta  Plasma,  Varenna, 
Italy,  Sept.  1977.   Also  ORNL  Report  ORNL/TM-6271,  1978. 


A  proof  of  existence  (not  convergence  of  the  numerical 
algorithm)  has  been  given  for  a  linearization  of  the  GDE  (6) 
[G.  Vigfusson,  Thesis,  NYU,  1977].   A  similar  linearization  is 
the  basis  of  the  entirely  2D  marching  algorithm  employed  by 
S.  C.  Jardin  (PPPL)  for  the  Grad-Hogan  model. 
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Six  Representative  lyD  Codes 


Code  Page  No, 


1.    Adiabatic,  2D,  Tokaraak,  Belt 

Pinch,  and  Doublet  (1974)   14 


2.    Adiabatic,  Axi symmetric 

Doublet  (1974)   35 


3.    Adiabatic  Tearing  -  2D  "Lens" 

(1975)   39 


4.    2D  Resistive  Tokamak  and 

Doublet  (1975-77)   41 


5.    Divertor  (1977)  59 


6.    Mirror  Machine  Code,  Simple  and 

Reversed  Field  (1978)  66 
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Code  1:   Adiabatic,  2D,  Tokamak, 
Belt  Pinch,  and   Doublet   (1974) 

I.    Purpose 

1.  To  test  adiabatic  theory  with  general  profiles 
for  Tokamak  compression. 

2.  To  test  generalized  adiabatic  theory  (with  change 
of  topology  and  jump  conditions  across  separatrix) 
for  Doublet. 

3.  To  test  numerical  efficacy  of  "lyD"  algorithm 
(formulated  1970)  for  solving  GDE's  by  alternating 
ID  determination  of  profiles  with  2D  evaluation  of 
geometric  coefficients  in  ID  equations. 

4.  To  test  theoretical  estimates  of  singular  current 
layer  induced  by  moving  separatrix. 

5.  To  study  general  effect  of  changing  shape  on 
pressure  and  current  profiles. 

6.  As  a  preliminary  to  coding  resistive,  lyD  Grad- 
Hogan  simulation  (where  similar  GDE's  occur). 
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II.   Description 


(a) 


(b) 


(c) 


1)  Scalar  pressure,  MHD,  arbitrary  6 
plasma  separated  from  conducting 
wall  by  vacuum  or  force  free  field. 

2)  2D  (with  third  component,  B  )  in 
rectangular  box  with  specified  BC 
on  li;  at  box  and  time  varying  coils 
between  plasma  and  box. 

3)  Pressure  profile  drops  to  zero  at 
plasma  vacuum  interface. 


4)  Simple  tokamak  or  belt-pinch  (a) ,  or  symmetric  Doublet  (b) , 
or  asymmetric  islands  (c) ,  including  adiabatic  transitions 
among  (a),  (b),  (c). 

5)  Plasma  evolution  governed  by  fixed  adiabatic  profiles 
(equivalent  to  entropy,  mass,  flux  and  rotational  transform, 
e.g.  "flux  conserving  Tokamak"). 

6)  Two  types  of  external  knobs  for  poloidal  field,  4)  (t)    at  wall 
and  time-varying  coils;  i>   represents  compression  and  coils 
control  shaping.   Also,  toroidal  vacuum  B  (t)  is  a  knob. 

Note:   Arbitrary  variation  of  external  poloidal  and  toroidal 

knobs  will  give  a  surface  current  at  plasma  edge;  since 
the  code  was  not  written  to  include  surface  currents, 
B   is  adjusted  to  follow  poloidal  field  to  maintain 
zero  surface  current. 
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III.   Equations 

A^  =  _  p  -  ff   ,   •  =  d/d^    ,       B^  =  f(ij^) 

'i^(x,y)  -^  ij^CV)  -^  V(iJ^)  ,  V  =  "volume"  inside  ^   =  const,  (area) 

p  =  u(>^)(i|j')^  defines  u,   '  =  8/9V 

f  =  v(i|j)ip'     defines  v 

li(4^)  and  v(ijj)   are  fixed  adiabatic  profiles  (p  and  f  vary 
with  constraints) . 

(1)  ^i)  =  -   \i{^^')^   -  yy  (4'')^"^'!^"  -  vv(ijj')^  -  v'^V  (GDE) 

(2)  <AiJj>  =  (KiJ;')'   ,   K  =  <|VV|^> 

(1)  is  the  complete  formulation  of  the  problem  (excluding  EC's); 
the  averaged  equation,  (2) ,  is  used  to  help  solve  (1) . 


Note:   by  setting  y   =    0,    the  code  will  construct  an  equilibrium 
with  specified  v(tj;)  and  p(4')  profiles. 
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IV.   Boundary  Conditions  and  Constraints 

1)  At  the  vacuum-plasma  interface,  ^,    Vtj;,  f  are  continuous 

[this  is  a  restriction  in  generality  -  cf.  (5)]. 

2)  In  case  (a)  ,  4^   at  V  =  0  (center)  and  ij^   at  plasma  edge 

P 

are  given  and  fixed  in  time  (i.e.,  for  varying  constraints); 
in  case  (c)  ,  i)      and  i)^   at  island  centers  are  fixed  as  well 
as  ij;  ;  in  case  (b)  '    ^i    -    ^2    """^  given.   In  an  adiabatic 
transition  from  (a)  to  (b)  or  (c)  ,  the  value  ^p      =   ^p^   =   \p 
is  preserved.   If  two  islands,  1  and  2,  are  shrinking, 
contours  in  regions  1  and  2  bearing  the  same  value  of  i^) 
reach  the  separatrix  simultaneously. 

3)  ij;  is  varied  for  compression  (vacuum  flux  =   "ip   -   ip      changes; 

P 

plasma  fluxes  ^      -   ip .    are  fixed). 

4)  I^  =  coil  strength  (fixed  locations  in  vacuum  regions) , 
varied  for  shaping. 

5)  The  solution  algorithm  described  below  determines  the 
plasma  shape  and  profiles  and  the  vacuum  poloidal  field, 
(|j(x,y);  the  vacuum  B   is  then  set  equal  to  f  at  the 
plasma  edge  (thereby  eliminating  any  surface  current) . 

6)  The  jump  conditions  across  Sep.  are  given  later. 
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V.    Basic  Algorithm 

To  solve  the  GDE,  (1),  requires  imposition  of  EC's  ^^    , 
ip,  ,  ^y   as  well  as  an  elliptic  BC  4^   (or  ijj)  .   Neither  the  left 
nor  right  side  of  (1)  dominates;  each  has  second  derivatives 
(in  2D  and  ID) .   The  conventional  elliptic  iteration  on  the 
RHS  will  not  work  -  for  one  reason  because  this  provides  no  way 
of  imposing  ^    .       The  average  GDE  (2)  is  a  second  order  ODE  for 
ij^  (V)  ,  assuming  that  K(V)  is  known, 

(3)   [K+yy  Cl^')^"'^  +  v^]ij;"  +  K'i|j'  +  iJ(ip')^  +  vv)(ijj')^  =  0 

This  PDE  can  take  EC's  at  each  end;  evidently,  in  case  (a), 

^{0)    =    ^         ,   ip(V  )  =  ^ 

The  distributed   inductance  K(V)  is  determined  by  the 
contours,  'P(x,y)  =  const,  [and  is  independent  of  the  profile, 
i>  (V)  ]  .   The  contours  are  found  by  inserting  an  approximation 
to  ^ {V)    on  the  right  side  of  (1)  and  solving  Poisson's  equation, 
It  is  advisable  to  eliminate  4^".  From  (3) 

^„  ^  _  K>'  +  y(4;')^  +  vvi|^'^ 
K  +  yui^') 

inserting  into  (1) , 


•)^-2   ^  v2 


(4)      A4;  =  F  S  K>'  [YU(4^')^"^  +  v^]  -  K[0(ij;')^  +  vv(<J^')^] 


The  form  (4)  is  used  in  2D  to  determine  the  contours. 
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The  iteration  starts  with  a  family  of  contours  and  K(V); 
(3)  solved  as  an  ODE  gives  i|^{V);  insertion  into  the  right  side 
of  (4)  gives  F(V);  the  previous  family  of  contours  converts 
this  to  F(x,y);  inversion  of  (4)  gives  a  new  set  of  contours. 

The  convergence  of  this  iteration  is  extremely  rapid. 
One  reason  is  that  the  contours  [hence  K(V)]  are  insensitive  to 
changes  in  profiles.   Another  reason  is  that  the  outer  loop  is 
almost  always  oscillatory,  allowing  critical  damping. 

For  an  initial  guess  at  the  contours  and  K(V)  it  has  thus 
far  been  found  adequate  to  take  Ai|i  =  1  in  the  entire  box. 
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VI.   The  Separatrix 

From  conservation  of  momentum  and  total  mass,  flux,  and 
energy,  the  following  relations  are  found 

i)   -rr   =  -rr  +    -rr  iv)   P   =  Pn  =  P-, 

^  f^      i)'       i)'  ^o   ^1   ^2 


ii)   v^  =  v^  +  v^  V)   f„  =  fi  -  f2 


iii)  y^^  -  y}/^  +  y^/^   vi)  K^^;  -  K^^'  +  K^^' 


\    =    'Pq/'I'I       '     (i  =  l,2)  ,  a^  +  a2  =  l 
(and  many  others) 

In  the  symmetric  case  Fig.  (b) ,  it  is  not  necessary  to 
be  aware  of  the  fact  that  islands  may  be  present  and  one  can 
still  use  the  simple,  one-region  algorithm.   Consider  an  inner 
region,  1+2.   For  a  given  (inner)  rp    contour  define  V{^)    = 
V^  +  V^  =  2V^   =  2V^.      Set  i>'     (inside)  "   J  ^[   =   J  ^2    (i-e.  use 
V  =  2V,  as  independent  variable).   V  and  i)'    are  now  continuous 
across  Sep.   Set  v(!|;)  (inside)  =  2v   =  2v  ;  y(^)  =  2^y,  =  ■^'^^7' 
y  and  V  are  now  continuous.   Moreover,  y(tjj)  and  v(tjj)  are  the 
same  invariant  functions  inside  and  outside;  they  are  not 
changed  by  a  moving  Sep.  and  are  known  a  priori  for  i|^  <    \p    <    ii    . 
This  procedure  loses  some  accuracy  by  glossing  over  the  singu- 
larity of  K(V)  at  V  =  V .   For  more  accuracy,  normalized 
coordinates,  e.g.  V/V    ,  are  used  in  each  region,  and  appropriate 
jump  conditions  connect  the  three  ODE's. 
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In  case  (c)  there  are  subcases  depending  on  whether  y  =   2 
or  not  and  depending  on  whether  conservation  of  energy  is 
utilized.   Not  all  cases  have  yet  been  run.   But,  in  the 
Doublet,  assuming  that  the  initial  state  is  simple,  and  the 
islands  are  created  later,  energy  is  conserved  and  all  relations 
(i)-(vi)  are  satisfied.  For  growing  islands  ("splitting"),  the 
originally  given  profiles  y   and  v   give  rise  to 

Pi  =  V^^       ,   i  =  1,2 
v.=va.   ,   i=l,2 

1      o  1 

For  shrinking  islands  ("mixing")  the  originally  given  profiles 
y .  and  v.  give  rise  to  y  and  v  as  in  (ii)  and  (iii) .  In  the 
splitting  case,  the  newly  created  profiles  are  not  independent 
of  path;  they  depend  on  a.  at  the  instant  of  splitting. 

The  basic  (one  region)  algorithm  is  changed  as  follows: 

1)  Given  the  contours,  we  have  three  profiles  K. (V)  and 
three  domain  volumes  V.(ZV.=V  ).  Three  ODE ' s  are  solved  subject 

J.       J-     JJ 

<l^„(v^)  =  4^i(v^)  =  4^2(^2)  =  -^s 

The  solutions  are  unique  if  4^   is  assigned.   The  value  of  i|) 
can  be  determined  in  several  ways,  e.g.  by  imposing  (i) ,  or 

Po  =  °'lPl  ^  "2P2 

(which  is  the  condition  that  energy  be  stationary) . 

2)  The  current  is  assigned  to  the  2D  mesh  separately  in 
each  region.   Then  the  Poisson  equation  is  solved. 
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In  an  iteration,  the  Geometry  step  (Poisson  solver)  will 
be  the  first  to  recognize  a  change  in  topology.   At  this 
signal  (or  slightly  later,  when  an  island  is  large  enough  to 
support  a  few  contours) ,  the  ODE  logic  is  changed. 
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VII.   Meshes  and  Interpolation 

1)  The  number  of  contours  in  a  simple  topology  or  in 
the  outer  region  (0)  is  fixed,  uniformly  spaced  in  i) ,   with 
the  option  of  subdividing  the  outermost  and  innermost  (center 
or  next  to  Sep.)  interval.   The  number  of  contours  in  an 
island  can  be  fixed  or  adjusted  by  a  formula,  according  to 
the  size  of  the  island. 

2)  The  ODE  mesh  is  in  this  early  code  based  on  the 
contour  mesh,  subdivided  by  a  factor  of  three;  there  is  no 
separate  interpolation  mesh. 

3)  Points  on  contours  are  located  by  successive  cross- 
ings of  mesh  lines  (discarded  if  too  close)  and  found  pre- 
cisely by  cubic  interpolation  of  four  colinear  meshpoints; 
linear  interpolation  is  very  poor  near  a  separatrix;  otherwise 
it  is  adequate. 

4)  The  previous  inversion  of  A(|j  =  F  is  preserved  as  a 
2D  array   4^£(x,y)  (tj^-label)  ;  4)^   is  more  convenient  than  V  to 
transfer  the  next  value  F  (V)  ->  F  {^ ^)    to  the  2D  mesh. 
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VIII.   Auxiliary  Algorithms 

1)   For  accuracy  (and  convergence)  when  the  plasma  is 
far  from  the  wall,  the  BC,  ip  =  ijj  ,  at  the  wall  is  carried  to 
the  plasma  edge  by  using  the  inductance,  L,  of  the  vacuum 
region  (between  the  plasma  and  the  box,  the  latter  indented 
by  the  coil  induced  separatrix.   It  is  not  necessary  to  draw 

contours  throughout  the  vacuum, 
only  one  vacuum  contour  immediately 
outside  the  plasma  (i|j^,V^,  etc.). 
Each  inversion  of  A(Jj  =  F  evaluates 
L^    for  the  current  geometry  as 


To  any  other  contour  inside  V.  the  inductance  is 


L  -  L^  + 


V 


dV/K   ,   (L  <  L^) 


V. 


The  plasma  edge  is  determined  by  choosing  V   such  that 

4>   -   i)      =  L  I 

p     p 

2)   The  2-point  boundary  ODE  is  solved  by  Picard  iteration 


using 


V  =  ij;' 


V  =  a  exp 


vdV  ^    \b      -  til 
^p    ^o 


•   •  Y— 2      • 
K+  yv'    +  W 


Y-2    2 

K  +  yyV   +  V 


dij; 


(to  determine  a) 
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IX.   Results 

1)  This  was  the  first  code  to  demonstrate  the  "Ip-D"  principle 
of  iterating  between  geometry  and  profiles. 

2)  This  was  the  first  code  to  compute  adiabatic  compression 

of  a  Tokamak  without  the  assumption  of  circular  contours 
[except  for  the  earlier  NYU  code  restricted  to  flat 
pressure,  simulating  ATC,  D.C.  Stevens,  Phys.Fl.  17_,  222  (1974)] 

3)  This  was  the  first  code  to  compute  the  evolution  of  a 
Flux  Conserving  Tokamak. 

4)  This  was  the  first  code  to  compute  the  time  evolution  of 
a  change  in  topology  (Belt  Pinch  to  Doublet) . 

5)  Compression  of  a  Belt  Pinch  was  found  to  circularize  the 
plasma  very  quickly  after  a  certain  critical  compression 
(as  had  been  shown  previously  for  a  flat  pressure  profile) . 

6)  A  sequence  was  found  in  which  the  plasma  was  observed  to 
first  bifurcate  (asymmetric  simple-contour  plasma  with 
symmetric  constraints),  then  islate. 

7)  This  code  demonstrated  the  theoretical  prediction  that 
the  rate  of  growth  of  islands  is,  under  certain  circum- 
stances, entirely  independent  of  dissipative  mechanisms 
(i.e.  "line-reconnection"  without  overt  dissipation). 

Results  are  reported  in: 

1)  P.N.  Hu,  H.  Grad,  D.C.  Stevens,  "Adiabatic  Compression  of 
MHD  Equilibria,"  Bull.  Am.  Phys.  Soc.  19^,  865  (1974). 

2)  D.C.  Stevens,  H.  Grad,  P.N.  Hu,  "Belt  Pinch  and  Doublet 
Compression",  Bull.  Am.  Phys.  Soc.  19,  865  (1974). 

3)  H.  Grad,  P.N.  Hu,  and  D.C.  Stevens,  "Adiabatic  Evolution  of 
Plasma  Equilibrium",  Proc.  Nat.  Acad.  Sci.,  72,  3789  (1975). 
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RUNS  -  Code  1  (19  74) 


A.  Adiabatic  Compression 

B.  Belt  Pinch  to  Doublet 

C.  Bifurcation 
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l.A-1 


Adiabatic  Compression  -  1 
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l.A-2 


Adiabatic   Compression    -    2 
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l.A-3 


Adiabatic  Compression  -  3 
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l.B-1 


Belt  Pinch   to   Doublet   -   1 
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l.B-2 


Belt   Pinch   to    Doublet   -    2 
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l.B-3 


Belt  Pinch   to   Doublet   -    3 
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B-4 


Belt   Pinch   to   Doublet   -    4 
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l.C 


Bifurcation 
Asymmetric  Equilibrium  with  Symmetric  Constraints 
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Code  2:  Adiabatic,  Axisymmetric  Doublet  (1974) 

This  code  is  an  axially  syminetric  version  of  Code  1.   It 
is  specialized  to  symmetric  islands  only.   It  was  written  in 
1974  by  D.C.  Stevens  and  M.S.  Chu  (GA) . 

I.    Purpose  -  see  Code  1. 


II.   Description 

1)  Scalar  pressure,  MHD,  arbitrary  3  plasma  separated 
from  conducting  wall  by  vacuum  field  and  external  coils. 

2)  Axially  symmetric  Tokamak,  Belt  Pinch,  or  Doublet 
with  arbitrary  cross-section,  aspect,  pressure  and  current 
profiles,  including  transition  from  Belt  Pinch  to  Doublet. 

3)  Fixed  adiabatic  profiles  with  pressure  zero  at 
vacuum  interface. 

4)  No  skin  current;  toroidal  field  is  adjusted  to  follow 
poloidal  field  variations. 


III.   Equations 

(1)      AiJ;  =  -  r^p  -  ff    ,    •  =  9/9^    ,    rB„  =  f  (4^) 


p  =  y(ijj)  {ip')^  ,         '  =  9/9V 

f  =  v(tJ^)ijj'/A      ,    A  =  <l/r^> 
GDE  as  in  Code  1  is  obtained  by  eliminating  p  and  f  in  favor  of 
y  and  v  in  (1) . 
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IV.  Boundary  Conditions  and  Constraints 

As  in  Code  1,  except  no  asymmetric  features. 

V.  Basic  Algorithms 

As  in  Code  1  except  that  both  K(V)  and  A(V)  are  taken 
from  equilibrium  contours  to  place  in  ODE  (averaged  GDE) . 

VI.  The  Separatrix 

Treated  as  described  in  Code  1  for  the  symmetric  case. 

VII.  Meshes  and  Interpolation 
As  in  Code  1. 

VII.  There  is  no  provision  for  handling  vacuum  region  by 
inductance. 

IX.   The  code  was  debugged  and  run  but  no  significant  parameter 
variations  were  made. 

Note  1:   The  main  conclusion  is  that  the  axially  symmetric  code 
is  written  and  runs  essentially  the  same  as  in  2D  for  symmetric 
islands.   The  asymmetric  code  (for  which  theory  predicts  quali- 
tative differences  from  2D)  was  also  coded,  but  it  was  not 
fully  tested. 

Note  2:   In  September  1975,  Code  1  was  again  converted  to  axial 
symmetry  during  a  visit  to  ORNL  by  D.  Stevens,  this  time  omit- 
ting the  Doublet  capability.   The  external  conducting  box  was 
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removed  (except  as  a  numerical  aid  to  satisfy  regularity  at 
infinity),  so  confinement  as  well  as  shaping  is  by  coils. 
Later  modifications  were  made  at  ORNL  by  J.  Hogan,  to  incor- 
porate dissipation  and  other  effects  and  by  D.  Nelson,  to 
investigate  dynamic  FCT  evolution. 

Note  3:   Recently  (1978) ,  W.  Grossmann  modified  Code  1  to 
include  the  axis  r=0  in  the  interior  of  the  plasma.   Moving 
the  outer  box  gives  a  simulation  of  liner  compression. 


ASPECT  RATIO  1.9  2 


ASPECT  RATIO  1.01 


AXIS 

Compression  of  Tokamak  to  Approximate  Spheromak 

(by  D.  Nelson) 
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Code  3:   Adiabatic  Tearing  -  2D  "Lens"  (1975) 

This  case  is  very  similar  in  general  conception  to  Code  1 
"Adiabatic,  2D,  Tokamak,  Belt  Pinch,  Doublet".   The  domain  is 
rectangular  but  this  time  periodic  in  x.   The  plasma  (finite  3) 


is  separated  from  the  walls  by  two  vacuum  regions.   All  equations, 
the  basic  algorithm,  and  the  boundary  and  jump   conditions  are 
the  same  as  Code  1.   There  are  two  major  subcases  for  the 
unperturbed  (slab)  equilibrium. 

a)  ^^    '    =  y      (finite  current) 

b)  \l)^    '    =    |y|    (sheet  current) 

There  are  three  major  differences  between  the  Doublet  and  Lens. 

1)  The  small  island  in  Code  1  (Doublet)  is  small  in  all 
dimensions;  in  the  Lens  it  is  small  in  only  one  dimension. 

2)  Growing  islands  in  the  Doublet  represent  splitting  of  region 
0  into  1  and  2;  here  a  growing  island  represents  mixing  of  1 
and  2  to  form  0  (and  vice  versa  for  shrinking  in  the  two 
cases) . 

3)  In  case  (a)  ,  the  adiabatic  profile  y  is  singular,  ]i   "^   l/ip 
for  small  ^. 

In  the  Lens  configuration  with  asymmetric  regions  1  and  2, 
conservation  of  energy  (i.e.  energy  stationary  with  regard  to 
perturbation  of  the  separatrix)  must  be  partially  discarded. 
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2 

Although  Vij;  is  continuous  and  p  +  1/2  f   is  continuous  across 

the  separatrix,  p  and  f  separately  are  not.   The  Lens  configura- 
tion, assumed  to  be  generated  by  a  perturbation  of  the  vacuum 
boundary  condition  is  "superstable"  with  regard  to  growing 
island;  the  potential  energy  would  look  like 


Growing  islands  (mixing)  follow  path  A.   Reversal  (shrinking  and 
splitting)  follows  path  B  which  is  reversible.   With  the  slopes 
and  curvatures  indicated,  there  is  an  adiabatic  energy  minimum 
at  equilibrium,  but  the  energy  is  stationary  on  only  one  side. 

Case  (b)  [with  nonsingular  y(i|^)]  offers  no  serious  dif- 

2 
ficulties.   Case  (a)  implies  \p   ^  V     near  the  center  of  the  island, 

This  means  B  is  very  small  near  the  center  and  the  field  topology 

is  extremely  sensitive  to  small  changes.   In  particular,  during 

iteration  the  topology  may  develop  strings  of  small  subsidiary 

islands  which  eventually  disappear  on  iteration.   There  are 

several  successful  but  intricate  methods  for  dealing  with  this 

problem.   Results  are  incomplete,  this  code  having  been  dropped 

to  proceed  with  resistive  codes. 

A  similar  axisymmetric  tearing  code  was  written  in  1975 

and  is  in  a  similar  state  of  suspension. 
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Code  4:   2D  Resistive  Tokamak  and  Doublet  (1975-77) 

I.    Purpose 

1.  To  construct  an  efficient  "1^"  transport  code  in 
the  Grad-Hogan  long  time  scale  formulation. 

2.  To  investigate  the  range  of  validity  of  Pfirsch- 
Schluter  diffusion. 

3.  To  investigate  the  coupling  between  geometry  and 
skin  effect,  geometry  and  plasma  diffusion,  skin 
effect  and  plasma  diffusion. 

4.  To  construct  an  efficient  MHD  code  with  real  (2D) 
geometry  as  a  foundation  for  next  generation 
"Kitchen  Sink"  Tokamak  and  Doublet  simulation. 

5.  To  compare  resistive  with  adiabatic  transition 
from  Doublet  to  Belt  Pinch. 

6.  To  investigate  current  profile  evolution  resulting 
from  geometry  and  topology  changes. 

7.  To  investigate  plasma  compression  on  the  skin  time 
scale . 

8.  As  a  preliminary  to  an  axially  symmetric,  finite  g 
transport  code. 
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II.   Description 


1.  Low  3,  scalar  pressure,  classical  resistivity,  MHD; 
temperature  profile  given  (no  energy  equation) • 

2.  2D,  plasma  separated  from  rectangular  box  by  vacuum  region 
with  coils. 

3.  Toroidal  and  poloidal  external  field  individually  varied. 

4.  Primarily  at  resistive  skin  time  scale  (slower  plasma 
diffusion  is  separately  computed) . 

5.  Tokamak,  Belt  Pinch,  Doiablet,  and  time  dependent  transitions, 

6.  Ability  to  handle  low  resistivity  (adiabatic  transition) 
with  accompanying  large,  localized  current  layers  (at 
plasma  edge  and  at  separatrix) . 


7.  Choice  of  plasma  edge  BC  with  skin  effect  and  skinless, 

8.  Specialized  to  y  =  2. 
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III.   Equations 

Low  3  formulation  [Nice,  19  70] 


||-      +    u-AiJ;    =    nAip    -    c(t) 

Aij;      =    J(iJ^,t)  [J   "    -    5^    '    P   "    P   +    I   ^^] 

<U.VUj>    =    U    ([»' 
o 

1         ^^ 
"o   =    ^(t)^      '         ^   =    -    i;     dt^ 

U   is  the  plasma  compression  induced  by  changing  the  zero- 
order  toroidal  vacuum  field,  f  (t) . 

r 
The  net  plasma  diffusion,  U,  =  o  u  'ds  {-    <u  •VV>)  is 

higher  order  on  the  skin  time  scale  and  can  be  computed 

(separately)  by  the  formulas 

V 


"gh^  Ss-^-^'^t  -  ^  J    (J^')^  dV  - 

V 

and 


V 

r 


V(JiJj')    dV         [Grad-Hogan] 


U        =   -   nKp'    -    cKijj'  [Pfirsch-Schliiter] 

The   complete    formulation   as    used    for   computation   is 

^^  +  u^tl^'   =  n  (KiJ^' )  '   -  c 

All;    =    J         ,  J   =     (Kip')  ' 
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IV.   Boundary  Conditions 

1)  A  poloidal  boundary  condition  ^    =   ijj(t)  [usually  constant  in 
space,  but  not  necessarily]  is  applied  on  the  external 
rectangular  box. 

2)  Coils  (in  this  code,  synunetrically  placed)  carry  varying 
currents  I  (t)  which  also  influence  the  poloidal  field. 

3)  The  toroidal  vacuum  field  f  (t)  can  be  varied  through  the 
convection  term  U  .   If  f   is  contant,  the  plasma  volume 

o         O  '  f 

V   is  constant  [otherwise  V  f   =  const. ] . 
p  p  o         ■■ 

4)  In  this  scaling,  poloidal  field  pressure  is  not  strong 
enough  to  compress  the  plasma;  changing  ijj  (or  I  )  will, 
however,  induce  a  skin  current  unless  a  special  skinless 
BC  (taken  from  the  adiabatic  problem)  is  used.   In  most 
cases,  one  of  the  following  BC's  has  been  used 

a)  \l)   -    given  at  box 

b)  Ki|j'  =  const.   at  V  =  V 

P 

c)  i)'    =    const.   at  V  =  V 

P 

BC  (b)  (constant  total  plasma  current)  gives  a  skin  current 
layer;  BC  (c)  can  be  shown  to  coordinate  (higher  order 
distinct  from  f  )  toroidal  and  poloidal  knobs  so  that  there 
is  no  skin  boundary  layer. 

5)  At  V  =  0  a  conventional  natural  BC  is  imposed  on  ^. 


6)    Treatment  of  Separatrix  is  described  separately. 
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V.    Basic  Algorithm 

For  a  given  K(V,t),  \p    is  marched  by  solving  a  conventional 
ID  diffusion  equation  in  0  <  V  <  V  with  conventional  BC  at 

V  =  0  and  choice  (usually  of  (b)  -  skin  or  (c)  -  skinless) 

at  V  . 
P 

In  an  early  version  of  this  code  [results  presented  at 
Madison,  1976]  ,  the  separatrix  was  ignored  and  the  two 
symmetric  islands  were  treated  as  described  in  Code  1,  taking 

V  =  V,  +  V„  =  2V, ,  as  independent  ID  variable  and  assuming 
that  ^   and  Kit'    are  continuous  across  the  Sep.  (essentially  by 
ignoring  its  location  or  existence) .   A  later,  more  accurate 
version  is  described  in  Sec.  VI. 

K(V,t.)  is  calculated  in  2D  at  intervals  t^  by  taking 
j(V,t.)  ^  J(x,y,t.)  from  the  diffusion  equation,  inverting 
Aij^j  =  J(x,y)  and  computing  a  new  K(V,t^)  from  the  contours  of 
^  I,  .      Between  t.  and  t .  .,  ,  K  is  assumed  to  be  linear  in  t. 
The  outer  loop,  assuming  that  all  quantities  are  known  at  t^, 
diffuses  to  t._|_-,  ,  recomputes  K^^,  ,  diffuses  again  from  t^  to 
1. 1  with  the  new  K  until  convergence. 

The  inner  loop  is  an  iteration  at  a  given  time  t^ ,  given 
J(V),  in  which  i) .    and  the  contours  are  iterated,  Ai/^  =  J[V(x,y)], 

J[V(x,y)]  ->-  i|)(x,y)      [Poisson] 

\1){x,y)      ->■  V(x,y)  [volume  within   contour] 
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VI.   Separatrix  Singularity  and  Normalized  Variables 

At  the  separatrix,  K  is  unbounded  as  log  | V-V  | ;  this 
is  nontrivial  since  K'  '^  1/ 1  V-V  |  occurs  in  the  inhomogeneous 
term,  J  =  {K\p')',    of  the  Poisson  equation.   The  actual  singu- 
larity can  be  evaluated  analytically  and  is  not  quite  so  bad. 
For  a  moving  singularity,  K  "^    log,  in  a  diffusion  equation 
"pj.  =  n(KiJ;')'  the  current  density  is  bounded  but  has  a  cusp 
almost  at  the  separatrix  (preceded  by  a  smooth  dip,  if  n  is 
sufficiently  small) .   The  height  of  the  cusp  (and  depth  of  dip) 
are  unbounded  as  n  becomes  small. 

The  following  normalized  independent  variables  are  taken, 

inside  and  outside  V 

sep 

V  =  V/V  0  <  V  <  1 


V  =  (V-V  )/(V  -V  )   ,     0  <  V  <  1 
s    p   s 

In  (V,t)  the  diffusion  equation  has  a  convection  term,  (V/V  )^' 
Both  V  (t)  and  K(V,t)  are  obtained  from  the  geometry  at  inter- 
vals t.,  K  by  linear  interpolation  in  t,  and  V  (t)  by  a  cubic 
spline  (so  that  V  (t)  is  continuous  at  t.)- 

A  moving  log  is  fitted  to  K  near  V  =  V  (t)  by  taking  a 

pair  of  contours  close  to  V  on  each  side.  K  =  a . +b .  log  I V-V  I 
'^  s  1   1      '    s ' 

is  fitted  to  each  side  (i  =  in,  out)  ;  b-  is  replaced  by  '^^i^i'^^n^ 
and  a.  refitted  to  the  last  contour  on  each  side.   The  log  fit 
does  not  make  much  difference  with  a  very  fine  2D  mesh,  but  it 
greatly  increases  accuracy  with  a  coarse  mesh  (which  cannot 
see  the  log) . 
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The  equilibrium  solver  is  very  sensitive  where  the 
current  peak  at  the  Sep.  is  very  large  [a  small  change  in  V 
makes  large  changes  in  J(x,y)].   The  outer  loop  is  found  to 
have  properties  similar  to  an  asymptotic  series  when  the 
current  peak  is  large  (small  r\   and  small  island)  ;  viz.  there 
is  an  optimum  outer  loop  time  step  t.  ,  -  t.;  if  too  large  or 
too  small,  many  iterations  are  required. 

In  a  transition  from  simple  geometry  to  islated,  the 
code  first  notices  a  change  in  topology  in  the  2D  Poisson 
inversion.   If  the  island  is  too  small,  it  is  ignored  (no 
contours  are  close  to  the  island)  and  K(0)  =  0  is  used  to- 
gether with  all  finite,  simple  contours. 

When,  after  convergence  of  the  geometry  at  t.,  the  island 
is  large  enough  to  fit  three  internal  contours  (plus  V  =  0  and 
the  Sep.)  K  is  computed  on  two  normalized  meshes  and  trans- 
ferred to  the  diffusion  mesh  before  diffusing  to  t .  ,  .   At 
the  same  time,  i)   is  recomputed  over  the  island  region  using 
(Kip')  '  =  J  and  the  new  K. 
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VII.   Meshes  and  Interpolation 

Two  separate  fine  meshes  are  used  inside  and  outside  the  Sep. 
for  ID  diffusion.   K,  evaluated  on  contours  (a  fixed  number 
outside,  varying  with  V   inside),  is  transferred  to  the 
diffusion  mesh  by  cubic  interpolation.   The  fine  diffusion 
mesh  is  also  used  for  transfer  of  J(V)  to  the  2D  mesh  by 
linear  interpolation.   The  contour  spacing  has  several  alter- 
natives: 

1)  specifying  the  number  of  contours  and  gap  from 
first  contour  to  Sep.  gives  a  geometric  spacing 
(constant  ratio  of  adjacent  gaps) ; 

2)  specifying  the  number  of  uniformly  spaced  contours; 

3)  optional  subdivision  of  "last"  subdivision  (near 

V  =  0  and  near  V  =  V  ) . 

P 

Although  possibly  unreasonable  at  first  glance,  the 
combination  of  ID  and  2D  meshes  allows  very  fine  resolution 
with  moderate  2D  mesh  and  contours  (say  32x64  and  10  outside 
contours)  in  which  a  very  sharp  current  profile  of  total 
width  less  than  one  2D  mesh  is  accurately  resolved,  and  an 
accurate  cusp  is  drawn. 

A  number  of  algorithms  for  further  improving  the  accuracy 
of  description  of  small  islands  have  not  yet  been  fully 
implemented . 
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VIII.   Auxiliary  Algorithms 

1)  The  outer  loop  is  almost  always  oscillatory  (this 
can  be  qualitatively  justified  theoretically) .   There 
is  therefore  a  great  advantage  in  convergence  to 
back  average  the  bare  result  of  iteration  with  the 
previous  iterate.   A  simple  exercise  in  control 
theory  adjusts  the  back  average  coefficient  in  terms 
of  the  observed  convergence  rate. 

2)  Back  averaging  is  used  in  both  outer  and  inner  loops; 
the  former  on  the  ID  array,  J(V),  and  the  latter  on 
the  2D  array,  ii  ^ .      The  second  also  serves  to  damp 
out  rapid  motions  of  the  separatrix. 

Not  yet  implemented: 

a)  Refining  2D  mesh  in  a  small  box  which  surrounds 
island  (to  allow  recognition  of  very  small  islands) . 

b)  Using  analytic  formulas  for  K(V)  in  small  island, 

extrapolating  back  to  precise  time  of  islation  with 

3/2 
analytic  growth  '^^  t    . 
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IX.   Results 


1)  Sample  results  including  a  demonstration  of  enormous 
deviations  from  Pf irsch-Schliiter  were  reported  at 
the  Sherwood  Theory  Meeting,  1975,  and  published  in 
the  Second  European  Conference  on  Computational 
Physics,  Garching,  1976. 

2)  Very  accurate  calculation  of  current  response  to 
deformation,  including  a  current  spike  precursor 
at  the  plasma  center  before  islation,  transient 
skin  effect  signature  at  the  center  of  the  plasma 

from  oscillatory  external  coils,  calculation  of 

2 

nJ  heating  at  plasma  center  from  oscillating 

separatrix  current,  and  accurate  reversible  (half 
sine  wave)  operation  at  small  n  were  reported  at 
Berchtesgaden  in  1976. 
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RUNS  -  Code  4:  (1975-77) 


A.  Transition  from  Belt  Pinch  to  Doublet 

B.  Pf irsch-Schluter  vs.  Grad-Hogan  Diffusion 
Velocity  Profiles. 

C.  Sharply  Peaked  Separatrix  Current 

D.  Growing  vs.  Shrinking  Island 

E.  Reversible  Run  with  Skin  Transient 
Residue  on  Axis. 
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(a)  Initial  State 


(c)   Fully  Developed  Islands 


4. A 


A:  Transition  from  Belt  Pinch  to  Doublet 


-53- 


\±\ 

_J 

CD 

ID 
O 
Q 


4.B 


-54- 


J23_£nPJLES_ 


Current   Profile 


rr.i57m  Irl.lOO 


SU3-  8.00 


4.C 


-55- 


4.D 
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External  coil  current  time  variation 


Tr  .0M266     1=  1.000  ^tW 


4.E-1 


E.l:   External  coil  current  at  maximum  with  plasma  current 
peak  at  separatrix. 
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E.2:   Separatrix  current  peak  partly  eaten  away. 


4.E-2 
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E.3:   AC  skin  transient  at  center  of  plasma 


4.E-3 
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Code  5:  Diver tor  (1977) 


This  is  a  modification  of  Code  4,  the  2D,  resistive,  low 
6  Tokamak.   The  geometry  is  as  shown.   There  are  four  coils 

in  the  vacuum,  the  left  pair 
stronger  than  the  right.   There 
is  a  simple  region  with  plasma 
(the  islands  are  in  vacuum) . 
The  edge  of  the  plasma  and  V 
are  determined  by  the  separatrix. 
^  =   V/V   serves  as  normalized  variable  just  as  V/V   in  Code  4. 
The  separatrix  represents  a  confluence  of  two  singularities, 
the  plasma  edge  and  the  separatrix.   The  analytic  solution 
near  the  X-point  has  a  log  and  is  extremely  difficult  to 
approximate  numerically  (e.g.  it  takes  a  2D  mesh  subdivision  of 
2     for  the  separatrix  angle  to  converge  from  72°  to  89°; 
90°  is  the  theoretical  prediction).   Nevertheless,  the  global 
plasma  behavior  can  be  accurately  represented  by  normal  means. 
The  essential  numerical  feature  of  this  code  as  compared  to 
Code  4  is  that  special  care  is  required  in  interpolating  near 
the  X-point. 
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Runs,  Code  5  (1977) 

The  following  cases  all  start  with  uniform  current  density 
in  the  plasma. 

In  Runs  A   to  D  the  ratio  of  currents  in  the  two  left 
corner  coils  to  the  right  coils  is  held  fixed  at  2.0;  the  four 
currents  rise  or  fall  in  proportion.   The  total  plasma  current 
is  held  constant  in  time  or  is  programmed  to  decrease  linearly 
in  time.   The  resistivity  is  either  constant  (n  =  -5  or  2)  or 
follows  the  profile: 

ri(v)  =    °-^^    ,   V  =  V/V  (t) 
1.1  -  V  ^ 

where  V   is  the  total  plasma  volume  (n  rises  by  a  factor  of  11 

ir 

at  the  edge) . 

A:  T]   -    2 ,    increasing  coil  current,  decreasing  plasma 
current  (two  skins  competing) 

B:  n (v)  profile,  increasing  coil  current,  constant 
plasma  current  (skin  competing  with  resistivity 
profile) 

C:   n  =  0.5,  constant  plasma  current,  reversible  coil 
current  (increases  to  a  maximum,  then  returns  to 
initial  value) 

D:   n  =  0.5,  constant  plasma  current,  increasing  coil 
current  (simple  skin  effect) 

These  results  were  presented  at  the  American  Physical  Society, 
Division  of  Plasma  Physics  Annual  Meeting  in  Atlanta,  November  1977, 
(Bull.Am.Phys.Soc. ,  Vol.  22,  p.  1095,  1977). 
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CoDE  5:   Initial  Configuration 


V 
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RuN  A  -  Final  State 


5. A 
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RuN  B  -  Final  State 


5,B 
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Code  6:   Mirror  Machine  Code, 
Simple  and  Reversed  Field  (1978) 

I.    Purpose 

1.  To  construct  an  efficient  transient  mirror  code  which 
can  realistically  take  into  account  all  long-time  scale  MHD 
effects  including  changes  of  topology,  and  also  allow  "Kitchen 
Sink"  physics  to  be  adjoined  as  in  Tokamak  simulation  codes 
(to  mimic  neoclassical  effects,  injection,  radiation,  etc.). 

2.  To  confirm  whether,  in  the  mirror  configuration,  one 
verifies  what  has  been  observed  in  other  configurations  (e.g. 
in  tearing  and  Doublet  simulations)  that  changes  in  topology 
and  island  formation  are  primarily  a  result  of  external  knob 
twiddling  (field  constraints)  and  are  only  secondarily  effected 
by  transport  or  stability.   In  particular  to  confirm  (what  was 
predicted  from  a  primitive  analytical  model  in  1974)  whether  a 
reversed  field  region  can  be  created  ideally  (i.e.  adiabatically) 
or  with  classical  transport,  solely  as  a  consequence  of  proper 
time  variation  of  the  external  mirror  coil  strength,  without 
injection.   The  code  is  constructed  to  be  neutral  with  regard 

to  topology  changes,  i.e.,  it  will  observe  whether  the  configura- 
tion is  simple  or  field-reversed  and  act  accordingly. 

3.  To  study  the  transient  physical  response  to  external 
field  coil  changes  and  injection,  to  study  the  long  term  decay 
or  approach  to  equilibrium,  to  observe  the  lifetime  of  reversed 
field  regions. 
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il.  Description 

1.  Scalar  pressure,  MHD,  arbitrary  3. 

2.  Axisymmetric,  periodic  in  z. 

3.  Plasma- vacuum  interface. 

4.  Simple  (a)  or  Reversed-Field  (b) ,  or  transition 
in  time  between  (a)  and  (b) . 

5.  Adiabatic  (generalized,  allowing  topology  change) 

or  classical  transport  -  scalar  resistivity  and 

2  2 

perpendicular  heat  conduction,  X  =  X/oa  x  . 

6.  External  knob:  time  varying  shaped  flux  BC  at  wall 
or  time  varying  mirror  coils  in  vacuum. 


»c^uu  r^ 


(a) 


(b) 
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III.      Equations 

A:      Grad-Hogan    (r,z,t) 


B   = 


J   = 


2-nr 

1 
2iTr 


2  1 

A*i|j   =    r   div(-^  Vi|j) 


1^  +   div(pu)    =    0 

^  +   u-Vp   +   YP   div  u  =    (Y-l)div(XVT)    +    (y-DnJ' 
9 1 

9i  +   u-vil;   =    nA*iJJ    -    c(t) 
t3t 

2    2 


n    =    n(T)     ,    X    =    X{T)/a)    T       ,    T   -    p/p 


p   =   pii))     ,    T   =    T(4>)     ,     p    =    P(t|^) 


A*,|j    =    -    4TT^r^(dp/d4j) 


B:      Averaged   Equations    in  V 


1^  <},{V,t)    =    *•     .    It    *^^'^^    "    *1 


f 
U   =   c)   u   dS    =    <u'VV> 
n 


p^   +    (pU) '    =    0 

P^   +    Up*    +    YPU'    -     (Y-1)(^*T')'    +     (Y-l)<nJ    > 

2         2 

i|j      +   Uijj'    =    4tt    n<r    >F    -    c 


F   =     (Kijj')  •  =-    p'/ip' 
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K  =  <|  Vv|^/4TT^r^> 


2       2    2   2 

<nJ  >  =  4tt  n<r  >F^ 


n  =  n-,_T 


-3/2 


X^   -  X^<r^>p^/T^^^{i>'}^ 


C:    Adiabatic  Variables 


M  = 


pdV   ,   0  <  M  <  1 


m  =  1  -  (1-M) 


1/2 


,   0  <  m  <  1 


1/2 
Pq  =  (l-M)-"/^  =  1  -  m 


^;;^  <|)(m,t)  =  (J)   ,   ^-j7  c!)(m,t)  =  ^ 


8m 


dt 


P  =  P/Pq 
*  =  2(^'/P 

^  =  4.^  +  U())'   ,  u  =  dv/dt 


a  ^  ±  ij;  =  4;Vp 


y  =  p/P 


Y 


m   is  adiabatic  independent  variable 


a  and  y  are  adiabatic  dependent  variables 
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D:    Complete  System  of  Equations 


(1-D)  ;    (1)   II  =  2TT2{n<r^>F) 


■^  ,  __  _  /s  . 


F  -  J  P(Kap) 


,2,   f^=A.B 


(Y-l)A. 


4P 


(Y-1) 


<r2>pV2 
22 

p  '     p  y 


B  =  4TT^(Y-l)n<r^>F^/p^ 


T  =  yp^'^'-'-Vp, 


(3)  i[Ko^   +   yyp^"^]  +P[a(Ka)*  +  yp^^"^^  -  0 


for  Y  =  2, 


(3) ■  p  = 


a(t) 


/Ka  +2y 


exp[-  2 


KP  +2y 


{2-D)  :    (4)  A*\i)   =  +  4Tr^r  F 


The  three  ID  profiles  are  a (m) ,  y (m) ,  p (m) ,  governed  by  (l)-(3) 
The  2D  geometry  is  obtained  from  (4) . 
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IV.   External  Boundary  Conditions 
At  wall,  r  =  a: 

il'  =  iJJ  [1  -  k  cos  z] 

^ it)    and  k(t)  are  given;  or  V  =  plasma  volume  is  given 
(instead  of  ip)  ;  or  ^=\l)      (k=0)  given  and  vacuum  coils,  I(t) 
(instead  of  k)  .   Normalize:  )|i  =  0  on  axis  in  case  (a)  ,    \p   =   0 
at  center  of  island  in  case  (b)  [to  correspond  to  adiabatic 
normalization  in  limit  A  ->  0  ,  n  -^  0]  . 

All  other  EC's  (at  edge  of  plasma,  at  r  =  0 ,  at  center 
of  island,  at  separatrix)  are  provided  internally  by  the 
structure  of  the  equations. 
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V.    Basic  Algorithm 

1)  a(m,t)  ,  y(m,t)  ,  p(m,t)  are  marched  in  ID  according 
to  two  coupled  diffusion-like  equations,  (1)  and  (2) ,  and  one 
ODE  (3) ,  used  to  carry  p  as  a  and  y  diffuse. 

EC's  are  required  for  (1)  and  (2)  at  m  =  0  and  m  =  1,  also 

for  one  adjustable  constant  in  (3). 

2 
During  this  diffusion  process,  <r  >  =  S(m,t)  and  K(m,t) 

are  assumed  to  be  given,  also  the  jump  in  Kijj'  across  the 

r 

separatrix  in  case  (b)  .   The  value  of  K;Jj'  =  o  Bdil  on  axis  or 
of  ijj  at  the  plasma  edge  is  transferred  to  the  2D  geometry 
calculation  from  diffusion. 

2)  At  intervals,  the  coefficients  S  and  K  are  computed 
from  a  2D  evaluation  of  ijj-contours;  also  the  jump  in  Ki|;'  across 
the  separatrix  in  case  (b)  is  evaluated  in  2D  and  transferred 
to  the  ID  diffusion. 

3)  Given  y (m)  ,  a(m)  ,  p (m)  after  diffusion,  the  function 
F(V)  for  computation  of  the  2D  equilibrium  is  obtained  from 

r 


V  =  2 


dm/p 


F  =  I  p(pKa) 


(F  is  obtained  by  differencing).   F(r,z)  is  placed  on  the  2D 
mesh  using  the  previous  V(r,z)  [actually,  a  tabulation  of 

'^£(r,z)  ]  . 

1 

4)   There  are  two  iteration  loops  (as  in  Code  4) . 
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5)  There  is  an  additional  fine  ID  "background"  or  inter- 
polation mesh.   This  is  to  separate  the  optimization  of  mesh 
size  for  diffusion  and  mesh  size  for  interpolation  and  transfer, 
All  transfers  between  ID  and  2D  are  made  via  the  background 
mesh. 

6)  At  present,  normalized  variables  are  not  used  in  the 
diffusion  equation,  but  they  are  used  in  transferring  data 
between  ID  and  2D. 

7)  There  is  provision  in  the  diffusion  algorithm  for  a 
variable  mesh-width  profile  (to  emphasize  boundary  layers) . 

8)  To  facilitate  code  changes,  the  diffusion  algorithm 
is  at  present  explicit  (Dufort-Frankel) . 
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VI.   Internal  Boundary  Conditions 

2 

In  the  simple  configuration,  <r  >  '\/  V  multiplies  each 

diffusive  second  derivative;  the  values  a(0)  and  p (0)  therefore 
float.   The  same  is  true  for  o   at  the  center  of  the  island; 

however  this  is  a  regular  point  for  y  which  therefore  requires 

2  2 
a  BC.   This  arises  from  X      =   X/w  t   where  co  'x.  B  '\>  0  at  the 

center.   Consideration  of  the  boundary  layer  introduced  by 

2  2  2  2 

replacing  w  t   by  e  +  w  t   gives  a  natural  boundary  condition 

relating  y  to  its  first  derivative  (this  is  equivalent  to 

A^T'  =  0) . 

At  the  edge  of  the  plasma,  a  lengthy  but  elementary 
analysis  leads  to  complicated  power  law  behavior.   Mathematical- 
ly, both  y  and  a  float  at  the  plasma  edge,  but  for  greater 
accuracy,  an  appropriate  exponent  can  be  fitted.   Heat  flow  and 
ohmic  heating  are  negligible  near  the  plasma  edge;  current 
density  can  be  small  or  large;  temperature  (in  the  absence  of 
radiation  cooling)  is  usually  large. 

At  the  separatrix,  for  the  adiabatic  problem,  K  and  Kijj ' 
are  bounded  and  nonzero  but  discontinuous.   In  the  resistive 
case,  1/ii'    and  K  are  logarithmically  infinite  and  moreover 
"discontinuous",  i.e.  the  coefficient  of  the  log  is  discontinuous, 

This  is  evident  from  the  dis- 
parate contours  inside  and  out- 

^1^1 — -^-i-  side. 

f  f 

l/i|;'  =  i  d£/B  ,  KiJ;'  =  o  Bd£ 

The  ODE  (3)  cannot  be  differenced  across  the  Sep. 
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It  can  also  be  seen  that  neither  the  a  or  y  diffusion  equation 
can  be  differenced  across  the  Sep.   Instead,  they  have  to  be 
matched.   The  jump  in  a  is  taken  from  geometry.   The  diffusion 

equation  for  a    is  treated  in  two  parts,  with  a  given,  moving 

2        .  .  . 

jump  in  a  and  with   <r  >F  continuous.   Similarly,  y  and  A^T 

are  continuous. 

Some  other  singular  behavior  at  the  Sep: 

2 

<r  >  'X'  1/log      [log  m  or  log  V  or  log  \\)] 

F  '^'  log  (current  is  unbounded) 

a  'V  i|j '  '^'  1/log 

At  present  these  weak  singularities  are  not  explicitly  noted; 
together  with  normalized  independent  variables  in  diffusion, 
they  will  be  made  explicit  in  the  future. 

An  interesting  remark  for  the  adiabatic  limit  is  that  ijj ' 


is  finite  and  continuous;  ©  dl/B   in  a  very  small  island  must 

r 
equal  cb  dil/B  on  the  entire  outer  axis! 
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Sample  Results,  Code  6  (1978) 

A:   In  a  long  run  with  the  mirror  coil  modulation  fixed 
one  expects  the  plasma  current  to  decay,  but  p  cannot  decay 
(in  the  absence  of  simulated  radiation  sinks  and  in  the  presence 
of  nJ  ) •   The  actual  behavior  is  found  to  be  an  approach  to  a 
flat  pressure  profile  with  a  surface  current  at  the  plasma  edge. 
Shown  are  the  flux  surfaces  and  current  profile  [precisely, 
F(V)  =  <J/r>]  for  the  times  early  and  late  in  the  run. 

B:   This  is  a  sample  calculation  with  resistivity  and  heat 
flow  turned  off  (adiabatic)  to  illustrate  the  formation  of  an 
island  in  a  simple  mirror  geometry.   The  flux  contours,  magnetic 
field  strength  on  axis,  and  current  profile  are  shown  at  the 
top,  for  the  beginning  of  the  computation.   As  the  mirror  ratio 
is  increased  (externally  applied  mirror  coil  currents)  the 
profiles  deform  as  shown.   As  in  earlier  doublet  calculations, 
there  is  a  current  peak  at  the  separatrix. 

References: 

1.  H.  Grad,  "The  Fat  Mirror  Machine",  Sherwood  Annual  Theory 
Meeting,  Berkeley,  Calif. ,  April  3-5,  1974. 

2.  H.  Grad,  W.  Grossmann,  D.C.  Stevens,  E.  Turkel,  S.  Wollman, 
"Macroscopic  Transport  Model  for  a  Reversed  Field  Mirror 
Machine",  Sherwood  Annual  Theory  Meeting,  Gatlinberg,  Tenn. , 
April  26-28,  1978. 

3.  D.C.  Barnes,  et  al. ,  "New  Results  in  High  Beta  MHD  Theory", 
paper  CN-L-3-1,  in  Proc.  7th  Intl  Conf.  on  Plasma  Physics 
and  Controlled  Nuclear  Fusion  Research,  Innsbruck,  Austria, 
Aug.  23-31,  1978. 
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late  in  run 


6. A 


A:   Current  sheet  develops  at  plasma  edge 
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This  report  was  prepared  as  an  account  of 
Government  sponsored  work.   Neither  the 
United  States,  nor  the  Administration, 
nor  any  person  acting  on  behalf  of  the 
Administration : 

A.  Makes  any  warranty  or  representation, 
express  or  implied,  with  respect  to  the 
accuracy,  completeness,  or  usefulness  of 
the  information  contained  in  this  report, 
or  that  the  use  of  any  information, 
apparatus,  method,  or  process  disclosed 
in  this  report  may  not  infringe  privately 
owned  rights;  or 

B.  Assumes  any  liabilities  with  respect  to 
the  use  of,  or  for  damages  resulting  from 
the  use  of  any  information,  apparatus, 
method,  or  process  disclosed  in  this 
report. 

As  used  in  the  above,  "person  acting  on  behalf 
of  the  Administration"  includes  any  employee 
or  contractor  of  the  Administration,  or 
employee  of  such  contractor,  to  the  extent 
that  such  employee  or  contractor  of  the 
Administration,  or  employee  of  such  contractor 
prepares,  disseminates,  or  provides  access  to, 
any  information  pursuant  to  his  employment  or 
contract  with  the  Administration,  or  his 
employment  with  such  contractor. 
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